Figures 1, S1-S3
Load packages
source("/Users/consti/Desktop/PhD/Statistical analysis/R-Scripts/Functions/summary.R")
require(maptools)
require(sp)
require(raster)
require(lattice)
require(ggplot2)
require(rgdal)
require(gridExtra)
Set directories and load data
#Environmental layers directory
raster.frost.dir = "/Users/consti/Desktop/PhD/Publication_material/24) Frost risk/R_data/CompleteRasters"
xtreme.years.dir = "/Users/consti/Desktop/PhD/Data/Global_layers/Environmental_layers/late_frost_maps/extreme_years"
#Continent shape file directory
continents.dir = "/Users/consti/Desktop/PhD/Data/Global_layers/Environmental_layers/continents"
oceans.dir = "/Users/consti/Desktop/PhD/Data/Global_layers/Environmental_layers/oceans"
land.dir = "/Users/consti/Desktop/PhD/Data/Global_layers/Environmental_layers/land"
#Output directory
output.dir = "/Users/consti/Desktop/PhD/Publication_material/24) Frost risk/R_output"
#-------------------------------------------------------------------------------------------------------------
#Create raster stacks and shape files
#####################################
#Late frost layers
complete.layers <- stack(paste(raster.frost.dir, "complete.rasters.grd", sep="/"))
complete.layers.dis <- disaggregate(complete.layers, fact=2, method="bilinear")
#Extreme frost years (only Northern Hemisphere)
xtreme.years.stack <- stack(list.files(path=paste(xtreme.years.dir), pattern='*.tif', full.names=TRUE))
xtreme.years.stack[xtreme.years.stack==0]=NA
xtreme.years.stack <- disaggregate(xtreme.years.stack, fact=4, method="bilinear")
#Land border shapefiles
continents <- readShapePoly(paste(continents.dir, "continent.shp", sep="/"))
oceans <- readShapePoly(paste(oceans.dir, "ne_50m_ocean.shp", sep="/"))
land <- readShapePoly(paste(land.dir, "ne_50m_land.shp", sep="/"))
Maps
Essential settings
# create own color schemes
col1 <- colorRampPalette(c("blue", "lightblue","yellow", "red", "red4"))
col2 <- colorRampPalette(c("red4","red4", "red3","red2","orangered","yellow", "lightblue","blue2", "blue3", "blue4","blue4"))
col3 <- colorRampPalette(c("tomato4", "tomato3","tomato2","tomato1","white", "skyblue1","skyblue2", "skyblue3", "skyblue4"))
col4 <- colorRampPalette(c("blue4", "blue3","lightblue", "yellow","orange"))
# continents and oceans polygons
pols1 <- list("sp.polygons", as(oceans, 'SpatialPolygons'), fill = gray(0.99), lwd = 0.1)
pols2 <- list("sp.polygons", as(land, 'SpatialPolygons'), fill = gray(0.8), lwd = 0.1)
1959-2016 average last frost date map
brks1 <- seq(0,190,0.5)
spplot(complete.layers.dis[["dateLastFrost0C"]],
at = round(brks1, digits=2),
maxpixels=2*ncell(complete.layers.dis[["dateLastFrost0C"]]),
col.regions = rev(col1(length(brks1))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"),
main=list(label='Average last frost date (<0C)',cex=1))
Figure S1a,b: 1959-2016 extreme LSF and extreme LSF change
#Maximum frost risk
DULF = complete.layers.dis[["Q95_DegreeZeroToFrostFour"]]
DULF[DULF>700] = 700 #set values above 700 to 700
brks1 = seq(0,700,0.5)
#Temporal trend (slope)
DULF2 = complete.layers.dis[["Slope_of_Global_DegreeZeroToFrostFour_Max_raster"]]*10
brks = seq(-50,50,0.1)
DULF2[DULF2> 50] = 50
DULF2[DULF2< -50] = -50
#Plot
figS1a = spplot(DULF,
at = round(brks1, digits=2),
maxpixels=2*ncell(DULF),
col.regions = (col1(length(brks1))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
figS1b = spplot(DULF2,
at = round(brks, digits=2),
maxpixels=2*ncell(DULF2),
col.regions = (rev(col3(length(brks)))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
pdf(paste(output.dir,"FigureS1ab.pdf",sep="/"), width=12, height=8)
grid.arrange(figS1a,figS1b)
dev.off()
grid.arrange(figS1a,figS1b)
Additional LSF maps (not shown in manuscript) and correlation with Fig. 1a
#get Pearson correlation coefficients
#####################################
#LSF 95% quantile vs LSF median
median_LSF = round(cor.test(as.vector(complete.layers[[8]]),as.vector(complete.layers[[4]]))$estimate,2)
#LSF degree-day threshol 0C vs 5C
dd5_LSF = round(cor.test(as.vector(complete.layers[[2]]),as.vector(complete.layers[[4]]))$estimate,2)
#LSF frost 0C vs frost -4C
extreme_LSF = round(cor.test(as.vector(complete.layers[[3]]),as.vector(complete.layers[[4]]))$estimate,2)
paste('Pearsons R =', median_LSF, "(median LSF),", dd5_LSF, "(dd5 LSF),", extreme_LSF, "(extreme LSF)")
## [1] "Pearsons R = 0.94 (median LSF), 0.94 (dd5 LSF), 0.86 (extreme LSF)"
##degreeDays 0 to frost 0, 50% quantile
DULF = complete.layers.dis[[8]]
DULF[DULF>650] = 650 #set value above 650 to 650
brks = seq(0,650,0.5)
#degreeDays 5 to frost 4
DULF1 = complete.layers.dis[[1]]
DULF1[DULF1>250] = 250
brks1 = seq(0,250,0.5)
#degreeDays 5 to frost 0
DULF2 = complete.layers.dis[[2]]
DULF2[DULF2>500] = 500
brks2 = seq(0,500,0.5)
#degreeDays 0 to frost 4
DULF3 = complete.layers.dis[[3]]
DULF3[DULF3>700] = 700
brks3 = seq(0,700,0.5)
#Plot
grid.arrange(
spplot(DULF,
at = round(brks, digits=2),
maxpixels=2*ncell(complete.layers.dis),
col.regions = (col1(length(brks))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"),
main=list(label='LSF, degree-day 0C, frost 0C, median',cex=0.9)),
spplot(DULF1,
at = round(brks1, digits=2),
maxpixels=2*ncell(complete.layers.dis),
col.regions = (col1(length(brks1))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"),
main=list(label='LSF, degree-day 5C, frost -4C, 95% quantile',cex=0.9)),
spplot(DULF2,
at = round(brks2, digits=2),
maxpixels=2*ncell(complete.layers.dis),
col.regions = (col1(length(brks2))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"),
main=list(label='LSF, degree-day 5C, frost 0C, 95% quantile',cex=0.9)),
spplot(DULF3,
at = round(brks3, digits=2),
maxpixels=2*ncell(complete.layers.dis),
col.regions = (col1(length(brks3))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"),
main=list(label='LSF, degree-day 0C, frost -4C, 95% quantile',cex=0.9)), ncol=2
)
Figure S3: Comparison between LSF change and extreme LSF change
#get Pearson correlation coefficients
#####################################
#LSF change (frost 0C) vs extreme LSF change (frost -4C)
extreme_LSF_change = round(cor.test(as.vector(complete.layers[[21]]),as.vector(complete.layers[[23]]))$estimate,2)
#LSF change (dd 0C) vs LSF change (dd 5C)
dd5_LSF_change = round(cor.test(as.vector(complete.layers[[19]]),as.vector(complete.layers[[23]]))$estimate,2)
#LSF change (95% quantile) vs LSF change (median)
median_LSF_change = round(cor.test(as.vector(complete.layers[[24]]),as.vector(complete.layers[[23]]))$estimate,2)
paste('Pearsons R =', median_LSF_change, "(median LSF change),", dd5_LSF_change, "(dd5 LSF change),",
extreme_LSF_change, "(extreme LSF change)")
## [1] "Pearsons R = 0.65 (median LSF change), 0.95 (dd5 LSF change), 0.12 (extreme LSF change)"
#Temporal trend 0C
DULF = complete.layers[[23]]*10
brks = seq(-50,50,0.1)
DULF[DULF> 50] = 50
DULF[DULF< -50] = -50
DULF[complete.layers[[15]]>=0.1] = 0
DULF = disaggregate(DULF, fact=2, method="bilinear")
#Temporal trend -4C
DULF2 = complete.layers[[21]]*10
brks = seq(-50,50,0.1)
DULF2[DULF2> 50] = 50
DULF2[DULF2< -50] = -50
DULF2[complete.layers[[13]]>=0.1] = 0
DULF2 = disaggregate(DULF2, fact=2, method="bilinear")
#Difference map
DULF3 <- complete.layers[[1]]
DULF3[DULF3<20000]=3
#both in right direction (at least one significant)
DULF3[complete.layers[[21]]>=0 & complete.layers[[23]]>=0 & (complete.layers[[43]]<0.1 | complete.layers[[41]]<0.1)] = 2
DULF3[complete.layers[[21]]<=0 & complete.layers[[23]]<=0 & (complete.layers[[43]]<0.1 | complete.layers[[41]]<0.1)] = 2
DULF3[complete.layers[[21]]>=0 & complete.layers[[23]]>=0 & complete.layers[[43]]<0.1 & complete.layers[[41]]<0.1] = 1
DULF3[complete.layers[[21]]<=0 & complete.layers[[23]]<=0 & complete.layers[[43]]<0.1 & complete.layers[[41]]<0.1] = 1
#both no significance
DULF3[complete.layers[[15]]>=0.1 & complete.layers[[13]]>=0.1] = 3
#both in wrong direction (at least one significant)
DULF3[complete.layers[[21]]>0 & complete.layers[[23]]<0 & (complete.layers[[15]]<0.1 | complete.layers[[13]]<0.1)] = 4
DULF3[complete.layers[[21]]<0 & complete.layers[[23]]>0 & (complete.layers[[15]]<0.1 | complete.layers[[13]]<0.1)] = 4
DULF3[complete.layers[[21]]>0 & complete.layers[[23]]<0 & complete.layers[[15]]<0.1 & complete.layers[[13]]<0.1] = 5
DULF3[complete.layers[[21]]<0 & complete.layers[[23]]>0 & complete.layers[[15]]<0.1 & complete.layers[[13]]<0.1] = 5
DULF4 = DULF3
#Plots
DULF3 = disaggregate(DULF3, fact=2, method="bilinear")
brks2 = seq(1,5,0.1)
figS3a = spplot(DULF,
at = round(brks, digits=2),
maxpixels = 2*ncell(DULF),
col.regions = (rev(col3(length(brks)))),
sp.layout = list(pols1, pols2),
colorkey = list(space = "right"))
figS3b = spplot(DULF2,
at = round(brks, digits=2),
maxpixels = 2*ncell(DULF2),
col.regions = (rev(col3(length(brks)))),
sp.layout = list(pols1, pols2),
colorkey = list(space = "right"))
figS3cUpper = spplot(DULF3,
at = round(brks2, digits=2),
maxpixels = 2*ncell(DULF3),
col.regions = (col4(length(brks2))),
sp.layout = list(pols1, pols2),
colorkey = list(space = "right"))
pdf(paste(output.dir,"FigureS3abc.pdf",sep="/"), width=9, height=11)
grid.arrange(figS3a,figS3b,figS3cUpper)
dev.off()
## quartz_off_screen
## 2
grid.arrange(figS3a,figS3b,figS3cUpper, ncol = 1)
# get percentages per continent and biome
DULFtable <- cbind(coordinates(DULF4), as.data.frame(DULF4), as.data.frame(complete.layers[[59]]))
DULFtable <- DULFtable[complete.cases(DULFtable$Q95_DegreeFiveToFrostFour),]
# continents intersection
coordinates(DULFtable) <- c("x", "y") #set coordinates
DULFtable$continent <- over(x=DULFtable, y=continents) # intersection
# converts the data to dataframe again
DULFtable <- as.data.frame(DULFtable)
DULFtable <- na.omit(DULFtable)
DULFtable <- DULFtable[!DULFtable$CONTINENT %in% c("Africa", "Antarctica", "Australia", "Oceania"),]
# get column with contradicting cases
DULFtable$wrong <- as.numeric(DULFtable$Q95_DegreeFiveToFrostFour %in% c(4,5))
DULFtable$continent_sub <- ifelse(DULFtable$CONTINENT=="North America", "North America",
ifelse(DULFtable$CONTINENT=="Europe", "Europe",
ifelse(DULFtable$CONTINENT=="Australia", "Australia",
ifelse(DULFtable$CONTINENT=="South America", "South America",
ifelse(DULFtable$x<80, "West Asia", "East Asia")))))
DULFtable$mix <- paste(DULFtable$continent_sub, DULFtable$WWF_Biomes_HalfDegree, sep="")
classification <- as.data.frame.matrix(table(DULFtable$mix, DULFtable$wrong))
classification$percent <- classification$`1` / (classification$`0` + classification$`1`)
classification$categories = factor(rownames(classification),
levels=c("Europe4", "West Asia4", "East Asia4", "South America4" ,"North America4",
"Europe5", "West Asia5", "East Asia5", "South America5" ,"North America5",
"Europe6", "West Asia6", "East Asia6", "South America6" ,"North America6",
"Europe8", "West Asia8", "East Asia8", "South America8" ,"North America8"), ordered=T)
classification$categories = factor(rownames(classification),
levels=c("Europe4", "Europe5", "Europe6","Europe8",
"West Asia4", "West Asia5","West Asia6","West Asia8",
"East Asia4", "East Asia5","East Asia6","East Asia8",
"South America4" ,"South America5" ,"South America6" ,"South America8",
"North America4","North America5", "North America6","North America8"), ordered=T)
figS3cLower = ggplot(classification, aes(x=categories, y=percent*100)) +
geom_bar(colour="black", stat="identity", position=position_dodge()) +
geom_hline(yintercept=50)+
coord_cartesian(ylim=c(5,95))+
labs(x = "", y = "Area different direction (%)") +
theme(legend.position= c(0.85, 0.8),
legend.background = element_rect(fill=NA,
size=0.5, linetype="solid"),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(angle = 90, hjust = 1),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = 'white', colour = 'black'))
pdf(paste(output.dir,"FigureS3c.pdf",sep="/"), width=11, height=4)
figS3cLower
dev.off()
## quartz_off_screen
## 2
figS3cLower
Figure S2: Extreme years
###################
#2007 North America
###################
DULF <- xtreme.years.stack[[16]]
DULF[DULF>300] = 300
DULF[DULF< -300] = -300
brks <- seq(-300,300,1)
#Change the spatial extent of stacks
boundary <- extent(-110, -70, 23, 45) #longitude min max, latitude min max
DULF <- crop(DULF, boundary)
spplot(DULF,
at = round(brks, digits=2),
maxpixels=2*ncell(xtreme.years.stack),
col.regions = (rev(col2(length(brks)))), # rev(terrain.colors(length(brks)))
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
###################
#2010 North America
###################
DULF2 <- xtreme.years.stack[[18]]
DULF2[DULF2>300] = 300
DULF2[DULF2< -300] = -300
brks <- seq(-300,300,1)
#Change the spatial extent of stacks
boundary <- extent(-90, -54, 39, 53) #longitude min max, latitude min max
DULF2 <- crop(DULF2, boundary)
spplot(DULF2,
at = round(brks, digits=2),
maxpixels=2*ncell(xtreme.years.stack),
col.regions = (rev(col2(length(brks)))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
############
#2011 Europe
############
DULF3 <- xtreme.years.stack[["Mean_Difference_of_DegreeZeroToFrostFour_2011_raster"]]
DULF3[DULF3>300] = 300
DULF3[DULF3< -300] = -300
brks <- seq(-300,300,1)
#Change the spatial extent of stacks
boundary <- extent(-5, 18, 38, 57) #longitude min max, latitude min max
DULF3 <- crop(DULF3, boundary)
spplot(DULF3,
at = round(brks, digits=2),
maxpixels=2*ncell(xtreme.years.stack),
col.regions = (rev(col2(length(brks)))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
############
#2013 Europe
############
DULF3 <- xtreme.years.stack[[20]]
DULF3[DULF3>300] = 300
DULF3[DULF3< -300] = -300
brks <- seq(-300,300,1)
#Change the spatial extent of stacks
boundary <- extent(2, 40, 44, 65) #longitude min max, latitude min max
DULF3 <- crop(DULF3, boundary)
spplot(DULF3,
at = round(brks, digits=2),
maxpixels=2*ncell(xtreme.years.stack),
col.regions = (rev(col2(length(brks)))),
sp.layout=list(pols1, pols2),
colorkey = list(space = "right"))
Temporal changes in frost probability
For regions in which frosts (< -4C or <0C) do not occur every year, these maps show the temporal change in the
probability of frost in a given year (red = increase in frost probability, blue = decrease in frost probability).
Change in frost probability is shown in percent per year.
frost_probability = complete.layers[[34:37]]*100
#crop upper and lower values
frost_probability[frost_probability > 4] = 4
frost_probability[frost_probability < -4] = -4
#Plot
brks <- seq(-5,5,0.1)
spplot(frost_probability[[base::c(2,4)]],
at = round(brks, digits=4),
maxpixels=2*ncell(frost_probability),
col.regions = rev(col3(length(brks))),
sp.layout= list(pols1, pols2),
colorkey = list(space = "right")
)
